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ABSTRACT 

By taking into account the local energy balance per unit volume between 
the viscous heating and the advective cooling plus the radiative cooling, we in- 
vestigate the vertical structure of radiation pressure-supported accretion disks in 
spherical coordinates. Our solutions show that the photosphere of the disk is close 
to the polar axis and therefore the disk seems to be extremely thick. However, 
the profile of density implies that most of the accreted matter exists in a moder- 
ate range around the equatorial plane. We show that the well-known polytropic 
relation between the pressure and the density is unsuitable for describing the ver- 
tical structure of radiation pressure-supported disks. More importantly, we find 
that the energy advection is significant even for slightly sub-Eddington accretion 
disks. We argue that the non-negligible advection may help to understand why 
the standard thin disk model is likely to be inaccurate above ~ 0.3 Eddington 
luminosity, which was found by some works on the black hole spin measurement. 
Furthermore, the solutions satisfy the Solberg-H0iland conditions, which indi- 
cates the disk to be convectively stable. In addition, we discuss the possible link 
between our disk model and ultraluminous X-ray sources. 

Subject headings: accretion, accretion disks — black hole physics — convection 
— hydrodynamics — instabilities 



Introduction 



The standard thin accretion disk model (IShakura fc Sunyaevl Il973[ ) has been widely 



applied to X-ray binaries and active galactic nuclei. Due to the basic assumption of the energy 
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balance between the viscous heating and the radiative cooling, such a model was known 
to be invalid for super-Eddington accre tion case, where the ad vective cooling is probably 
significant. Instead, the slim disk model (lAbramowicz et al.lll988l ) was introduced to describe 
super-Eddington accretion disks. However, there exists some conflict between the theory 
and t he observation. The theory predic ts that the advection is negligible for L < Lead 
Watarai et all l2000l ; ISadowskil 120111 ) , where Z/Edd is the Eddington luminosity, which 



(e.g. 

indicates that the standard disk model should be valid up to L Edd . On the contrary, some 
works on the black hole spin measure ment showed that the s tandard disk model is likely 
to be inaccurate for L > 0.3LEdd (e.g., 



McClintock et a 



general model for optically thick disks (e.g., ISadowski 



2011 



2006) . Moreover, even th e recent 



Sadowski et al.ll201ll ). which 



unifies the standard thin disk and th e slim disk, could not help to obtain a self-consistent 
spin parameter for L > 0.3L E dd (e.g.. IStraub et al.ll201ll ). In our opinion, the above conflict 



may be resolved if the vertical structure is well incorporated. 

Most previous works on accretion disks focused on the radial structure in cylindrical 
coordinates (R, <p, z). For the vertical structure, however, a simple well-known relationship 
U H = Cg/fix" or "iJ^x/cs = constant" was widely adopted, where H is the half- height of 
the disk, c s is the sound speed, and Qk is the Keplerian angular velocity. Such a relationship 
comes from the vertical hydrostatic equilibrium with two additional assumptions. One is 
the approximation of gravitational potential: i()(R, z) ~ ip(R, 0) + f2|-z 2 /2, and the other is 
a one - zone appro ximation or a polytropic relation p tot = K,p l+l l N in the vertical direction 
(e.g.. lHoshilll977l ). where p tot is the total (gas plus radiation) pressure and p is the density. 



Obviously, the above assumptions work well for geometrically thin disks, but may be inaccu- 
rate for the mass accretion rate M approaching the Eddington one MEdd, for which the disk 
is probably not thin. Consequently, the relationship u HQk/c s = constant" may be invalid 
for M > M Edd . 

Without the potential approximation, our two previous works investigated the geomet- 
rical thickness o f accretion disks and the validity of the relationship u HQk/c s = constant". 
Gu fc Lul (120071 ) adopted the explicit gravitational potential in cylindrical coordinates and 
found that the abov e relationship is inaccurate for M > M Edd , and therefore the disk can be 
geometrically thick. iGu et al.l ( 120091 ) took spherical coordinates to avoid the approximation 
of gravitational potential, and found that an advection-dominated accretion disk is likely 
to be quite thick. In these two works, however, the polytropic relation is still adopted in 
the vertical direction, which takes the place of the energy balance per unit volume between 
the viscous heating and the advective cooling plus the radiative cooling. The validity of 
such a polytropic relation, however, remains questionable, in particular for large M due to 
dominant radiation pressure. 
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The purpose of this paper is to revisit the vertical structure of radiation pressure- 
supported disks by taking into account the local energy balance and to study the variation 
of energy advection with mass accretion rates. The paper is organized as follows. Equations 
and boundary conditions are derived in Section 2. A global view of the solutions in the rh-r 
diagram is presented in Section 3. For a typical radius r = 10r g , the vertical structure and 
the energy advection are investigated in Section 4. The two-dimensional solutions and the 
convective stability are studied in Section 5. Summary and discussion are made in Section 6. 



2. Equations and boundary conditions 
2.1. Equations 



We consider a steady state axisymmetric accretion disk in spherical coordinates (r, 9, <f>) 
and use the Newtoni an potential, if) = —GM/r, where M is the black hole mass. Following 
Narayan fc Yil (119951 ). we assume vg = for simplicity whi ch means a hydrostatic equilibrium 
in the 9 direction. Simulations, (e.g.. lOhsuga et al.ll2005l Figure 3), however, revealed that 
vq will be significant for extremely high accretion rates such as M = 1000LEdd/c 2 . As shown 
in the following sections, our solutions mainly correspond to M around MEdd- F° r such 
accretion rates, the validity of v g — remains a question. 



The basic equations of continuity and momentum take the forms (e.g.. lKato et al.ll2008l ): 
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where v r and are respectively the radial and azimuthal velocity, F r and Fg are respectively 
the radial and vertical radiation flux, p is the gas pressure, K es is the opacity of electron 
scattering, and T rt j, is the r<p component of the viscous stress tensor, t t $ = uprd(v^/r)/dr. 
Following the spirit of a stress prescription, we assume the kinematic viscosity coefficient 
v = aclr/v K , where c s is the sound speed defined below (Equation (7)), and v K = {GM/r) 1 / 2 
is the Keplerian velocity. 



We would stress that, even though the a stress prescription has been widely adopted for 
theoretical studies, simulations of magnetorotational turbulence have shown that the stress 
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does not well scale locally with the pressure. For instance, the simulations on thin disks by 
a shearing box showed that, the time and box-averaged results are likely to support that the 
stress is proportional to the therma l (gas plus radiation) pressure (e.g., iHirose et al.ll2009al 
Figure 3). However, Figure 16 of IHirose et al.l ( l2009bl ) shows that the maximal thermal 
pressure is located on the equatorial plane, whereas Figure 11 shows that the maximal stress 
is obviously not at the same place. These two figures reveal that the stress is not proportional 
to the pressure locally. In the present study, for simplicity, we will keep the local a stress 
prescription for numerical calculation, which is a weak point of this work. 



The energy equation including gas and radiation is written as (e.g.. lOhsuga et al.ll2005[ ) 

V ■ [(e + E)v] = -pV ■ v - Vv.P - V ■ F + $ vis , (5) 

where e and E are the internal energy density of the gas and the radiation, respectively. 
P = fE is the radiation pressure tensor, $ vis is the viscous dissipative function, and the 
radiation flux F is expressed as 



Xc 



-VE 



(6) 



In this work, we focus on the region inside the photosphere, so we can take the well-known 
Eddington approximation, i.e., A = 1/3 and the Eddington tensor / = 1/3. 

Since we only study the radiation pressure-supported disks, the gas pressure p and the 
gas internal energy density e will be dropped in our calculation. In order to avoid directly 
solving the partial differential equations, some assumptions on the radial derivatives (d/dr) 



are required. Follow ing the spirit of self-similar assumptions (e.g., iBegelman fc Meierlll982 



Narayan fc Yilll995[ ). we adopt the following radial derivatives for c s and E: 



<91nc s 
<91nr 

where the sound speed c s is defined as 
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Based on the above two radial derivatives, the following four derivatives can be inferred from 
Equations (l)-(7): 
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With all the above derivatives, we can remove the "<9/<9r" terms in Equations (2) and 
(4)-(6), and the following equations are then obtained from Equations (2)-(6): 

1 „ 5 



2,^2, 2 
2«r + + 



J K 
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v}c0te = -^F e , (9) 
= (10) 

-^J-^)=S?I(^W). (11) 
f = -^f 5 ^ • < 12 > 

* = }|— ■ (13) 
6 rp/t cs 

The seven equations, Equations (7)-(13), enable us to solve for the seven variables: v r , 
V0, c s , p, E, F r , and Fq. There are two differential equations in this system. In addition, the 
position of the surface is unknown. Thus, totally three boundary conditions are required to 
determine a unique solution. 



2.2. Boundary conditions 

An obvious boundary condition on the equatorial plane is F e = 0. However, this condi- 
tion is not applicable for numerical calculation since it is automatically matched as indicated 
by Equation (9). Combining Equations (9) and (11) we can derive the following equation: 

cote±(vl) = vl+ r -^(3pvl-E). 

An alternative boundary condition on the equatorial plane is then obtained from the above 
equation (the left-hand side is zero thus the right-hand side should also be zero): 

vl + r -^(3pvl-E) = (0 = |). (14) 

The second boundary condition is a definition of the surface. We define the photosphere as 
the position above which the optical depth is around unity. The condition can be written as 

res = TK cs p 2 (^) 1 = 1 (0 = Bo) , (15) 

where #o (0 < #o < tt/2) is the polar angle of the photosphere. The third condition is related 
to the mass accretion rate: 

M = -2irr 2 / pv r sm6 d6 . (16) 
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3. Solutions in m-r diagram 

In our calculation we set M = lO Mm, k rs = 0-34 cm 2 ^" 1 , and a = 0.02, where the value 



of a is taken from recent simulations ( IHirose et al.ll2009a| ). The Eddington accretion rate is 



expressed as M Edd = 47iGM/r]CK es , where r) is a radiative efficiency of the flow. We choose 
1] = 1/16 since it is comparable to the Schwarzschild black hole efficiency of 0.057. The 
dimensionless accretion rate is defined as m = M / M-&aa- 

With the equations and boundary conditions in Section 2, we can numerically derive 
the ^-direction distribution of physical quantities for a given m at a certain radius r. The 
radiation pressure-supported disk solutions in the m-r diagram are shown in Figure 1, where 
r g = 2GM/ c 2 is the Schwarzschild radius. The parameter space is divided into three regions 
by two parallel solid lines, roughly with m oc r. The region above the upper solid line is 
denoted by "Outflow" , where we cannot find solutions. No solution exists probably due to 
the assumption of vq = in advance. In our view, a real flow located in this region may have 
v e 7^ and the inflow accretion rate may decrease inward. The physical understanding could 
be that, for high accretion rates and particularly for the inner radii, the viscous dissipation 
may be sufficiently large such that the radiation pressure is too strong to be balanced by 
the gravitational force. Thus, outflows may be driven by the radiation pressure and the 
inflow rh drops inw ard. On the other hand, simulations of supercritical accretion flows (e.g., 



Ohsuga et al.ll2005l Figure 6) showed that the inflow accretion rate roughly follows the m oc r 
relationship for m = 1000LEdd/c 2 at r out = 500r g (corresponding to m = 62.5 due to the 
definition of M-^dd with rj = 1/16). The slope of the upper solid line in Figure 1, which may 
be regarded as maximal accretion rates due to our calculation, agrees well with the slope in 
the above simulations. 

The region under the lower solid line is denoted by "Gas pressure" , where no solution is 
found either. In our understanding, it is probably because the gas pressure cannot be ignored 
in this region, which may be in conflict with the radiation pressure-supported assumption. 
We would point out that, the lower solid line in this diagram is higher than the well-known 



line w hich separates the inner and middle regions of standard thin disks (IShakura fc Sunyaev 



19731 ). The reason is that, the gas and radiation pressure are comparable for the latter, 
whereas the radiation pressure-supported disk may require the accretion rate to be higher 
such that the radiation pressure sufficiently dominate over the gas pressure, and therefore 
the effect of gas pressure on the vertical structure can be completely ignored. 

The region between the two solid line, denoted by "Radiation pressure" , which means 
that the radiation pressure is completely dominated, corresponds to the solutions of our 
main interest in this work. In Section 4, we will focus on the vertical structure and the 
energy advection at a typical radius, r = 10r g , as indicated by the vertical dashed line in 
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Figure 1. In Section 5, we will study the two-dimensional solutions for a typical accretion 
rate rh = 0.6 in the range 6r g ^ r ^ 12r g and < 9 ^ 7r/2, as indicated by the horizontal 
dot-dashed line. In addition, we would point out that for inner radii such as 3 ~ 5r g , the 
two solid lines in Figure 1 may deviate from a real black hole accretion system due to the 
Newtonian potential used in this work. 



4. Solutions at a typical radius r = 10r g 

4.1. Vertical structure 

In this section, we will focus on the solutions at a typical radius r = 10r g . Figure 2 shows 
the vertical structure of the disk with m = 0.6. In Figure 2(a), the dot-dashed, dotted, solid, 
and dashed lines show the vertical distribution of the dimensionless density (p/p ), radial 
velocity (v r /vjt), azimuthal velocity (v^/vg), and sound speed (c s /vk), respectively, where po 
is the density on the equatorial plane. It is seen that p significantly decreases, whereas c s and 
| i; r | increases, from the equatorial plane to the surface. In Figure 2(b), the solid line shows 
the variation of r es (defined in Equation (15)), where the photosphere (r es = 1) is located at 
9 m 4°, quite close to the polar axis. The disk seems to be extremely thick according to 
the position of the photosphere. However, the profile of p implies that most of the accreted 
matter exists in a moderate range around the equatorial plane, such as 7r/4 < 9 < 3% /4, 
which is more clear in Figures 4 and 5 (discussed below). The dashed line shows the variation 
of \F e \/cE. There exists \F e \/cE < 1/3 for the whole solution, which indicates that the 
Eddington approximation is valid and the solution is therefore self-consistent. 



With a more general viscosity, iBegelman fc Meierl (119821 ) studied a geometrically thick 



radiation pressure-supported model for supercritical accretion disks. They showed that there 
exists a narrow empty funnel along the rotation axis with a half-opening angle < 4°. 6. As 
seen in our Figure 2(a), the density drops sharply close to the photosphere, thus a nearly 
empty funnel also s eems to exist in our model. The difference is that, the disk surface in 



Begelman &: Meierl (119821 ) is the position where some physical quantities such as v r diverges, 
whereas our surface is defined as the position where Equation (15) is matched, and no 
divergence appears in our solutions. 

For a real disk with rh = 0.6, the photosphere may exist between 9 = 45° and the present 
result (« 4°). Our argument is as follows. There are two possible reasons that may cause the 
present photosphere quite close to the polar axis. One is that we have ignored the radiation 
force from one side (e.g., 9 = 9 and < < n) to the other (e.g., 9 = 9 and 7r < < 2%). 
The other reason is that we consider only the rcf) component of the stress tensor, which may 
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cause inaccurate results for small such as strong shearing of the angular velocity f2 in the 
vertical direction, where Q = v<p/(rsm6). Nevertheless, the r<p component assumption may 
work well for moderate 9: tt/A < 9 < 3n/4. The fact that the surface condition could not 
be matched in the range n/A < 9 < 3n/4 indicates that the half-opening angle of the disk 
(tt/2 — 9q) is likely to be larger than 7r/4. 

The profile of c s in Figure 2(a) is quite differe n t from that in the previous works with 



a vertical polytropic assumption (e.g., iHoshil 119771 ; iGu et al.l 120091 ). Under the polytropic 
relation p tot = JCp 1+1 / N (normally 1.5 < N < 3, and for radiation pressure-dominated case, 
Ptot can be replaced by E/3), c s will decrease continuously from the equatorial plane to the 
surface. The reason for the opposite behavior of c s , as implied in Figure 2(a), is that p drops 
faster than E from the equatorial plane to the surface. Figure 3 shows the variation of the 
quantity din E /din p with 9 for m = 0.5 (dashed line), fn = 0.6 (solid line), and m — 1 
(dotted line). If the polytropic relation works well, din E/ din p should be a constant of 
1 + 1/N. It is clearly shown in Figure 3 that, however, d In E/d In p varies significantly with 9 
rather than being a constant. More importantly, din E/ din p < 1 indicates that N is negative 
thus unacceptable. We therefore argue that the polytropic relation should be unsuitable for 
describing the vertical structure of radiation pressure-suppo rted disks. Moreover, sin ce the 
energy advection is relevant to c s (e.g., Q a dv — Mcl/2nR 2 in lAbramowicz et al.lll995l . where 
Qadv is the advective cooling rate per unit area), we may expect essentially different results 
on the strength of advection. 



4.2. Energy advection 

Figure 4 shows the variation of the vertically averaged advection factor / adv with the 
mass accretion rate m, where / a d v is defined as / a d v = Q&dv/Qvis- The quantities Q a d v and 
Q v i S are expressed as follows: 

Qadv = r q adv sin 9 d9 , (17) 

Je 

/•TT—00 

Qvis = r q vis sin 9 d9 , (18) 

Je 

where g a d v = —v r E/2r and g v i S = —3pv r v^/2r are respectively the advective cooling rate and 
the viscous heating rate per unit volume, as implied by the left-hand side of Equation (11). 

The solid line in Figure 4 corresponds to the total accretion rate integrating from 9$ 
to 7r — 9 , as shown by Equation (16), whereas the dashed line corresponds to the specific 
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accretion rate integrating from 9 = 7r/4 to 9 = 3tt/A, i.e., 

/■37T/4 

= -2vrr 2 / /w r sin 9 d9 . 

Jtt/4 



(19) 



The reason why we calculate for M^u is that the rip stress assumption may work well for 
7r/4 < 9 < 37r/4. As shown by the horizontal range of the solid and dashed lines, most of the 
accreted matter exists in this specific range, e.g., m^u = 0.52 corresponding torn = 0.6. The 
figure also shows that / adv rapidly increases with increasing m in the range 0.5 < rh < 1.1. 
More importantly, the value of / a d v (0.2 < / adv < 0.8) indicates that the energy advection is 
significant even for sub-Eddington accretion disks. 

Such a result is quite different from the previous one, where advection was found to 
be significant only for super-Eddington accretion case. IWatarai et al.l (120001 ) introduced an 
elegant formula to describe the M — L rela tionship based on their nu merical solutions under 
the well-known Paczyhski-Wiita potential ( jPaczyriski fc Wiitalll980l ). Their Equations (15)- 
(19) imply that, for the position r = 10r g , advection is negligible for M < 67L Edd /c 2 
(~ 4M Edd ). For the whole disk, advection is negligible for M < 20L Edd /c 2 (1.25M E dd)- 
Such a critical M for the whole disk was confirmed by some recent global solutions under 
the general relativity. Figure 4.11 of ISadowskil ( 1201 ll ) shows that advection is negligible 
for L < I/ Edd for any spin parameter a*. For a* = 0, i.e., the Schwarzschild black hole, the 
critical M is just around M Edd . In our opinion, the different results on the advection between 
the above two works and ours are related to the different approach to describing the vertical 
structure. 



A significant difference is that IWatarai et al.l (120001 ) and ISadowski et al.l ( 1201 ll ) chose 
the cylindrical coordinates whereas we adopt the spherical coordinates. Of course, the 
final results should not depend on the coordinates used. However, as pointed out by 
Abramowicz et al.l ( 119971 ). there are some interesting differences between the equations writ- 
ten in cylindrical and spherical coordinates. There is no centrifugal force in the z direction in 
cylindrical co ordinates, whereas t here is no gravitational force in the 9 direction in spherical 
coordinates. lAbramowicz et al.l ( 119971 ) claimed that it is exactly this property that makes 
the spherical coordinates much better adapted for describing the flow near the black hole 
horizon. Here, we argue that the spherical coordinates should be more suitable for describing 
geometrically thick disks as follows. In cylindrical coordinates , in the z direction, whether 
with a polytropic relatio n between the pressure and the density (IWatarai et al.ll2000l ). or with 
the local energy balance (ISadowski et al.ll201ll ). an approximation for the gravitational force, 
i.e., dip/dz = &k z i was adopted for describing the vertical structure. Such an approximation 
will probably be invalid for z/r > 1. In particular for z — > oo, the approximate force goes 
to infinity whereas the real force ought to vanish. Thus, the cylindrical coordinates seem 
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unsuitable for studying geometrically thick disks. In other wo rds, a geometrically th ick disk 
solution in cylindrical coordinates may not be self-consistent. ISadowski et al.l ( 120111 ) limited 
their solutions by M ^ 2M Edd probably due to this reason. As shown by their Figure 10, the 
maximal value of H/r for M = 2Msdd is ~ 0.4. For higher M, the value of H/r will be even 
larger thus the solution based on the approximate force may be inaccurate. On the contrary, 
in spherical coordinates, there is no need to make approximation for the gravitational force. 
The centrifugal force in the 9 direction, which takes the place of the z-direction gravitational 
force in cylindrical coordinates, is derived in this work by solving the vertical differential 
equations. Thus, our approach to the vertical structure seems to be more reasonable. 



We would agree that the solutions in lSadowski et al.l (120111 ) are likely to be self-consistent 
since their H/r is significantly less than unity, in particular for the solutions with M < MEdd- 
Then what are the reasons for the quantitative difference in the advection for M < Mgdd 
between their solutions and ours? In our under standing, there exist th ree possible reasons 
as follows. First, as mentioned in Section 3 of ISadowski et al.l (120111 ) for their numerical 
methods, the vertical structure is derived by a given advection factor / a d v in advance. The 
value of / a dv is probably obtained by solving the radial structure on the equatorial plane. 
Moreover, their / a d v is assumed to be uniform in the z direction. On the contrary, we obtain 
a varying / adv by solving the vertical equations. As shown by our Figure 8, / a d v increases 
significantly with z. We can therefore expect that the vertically averaged / a d v at a cylindri- 
cal radius will also be sig nificantly larger than that at z = 0, which may expla in why our 
/ a d v is larger than that in lSadowski et al.l (120111 ). Second, ISadowski et al.l (120111 ) assumed a 
uniform vr and in the z direction, and used the Keplerian strain to calculate the viscous 
dissipation. In our method, however, we include varying v r and in the ve rtical direction, 



and th e viscous dissipation is calculated based on v<f, instead of vk- Third. ISadowski et al. 
(120111 ) assumed v z = whereas we have vg = 0. As stressed by lAbramowicz et al.l (119971 ). 
since the stationary accretion flows resemble quasi-spherical flows (9q constant) much more 
than quasi- horizontal flows (H constant), vg = may be a more reasonable approximation 
than v z — 0. Moreover, the above three reasons may also be responsible for the different 
results in the convective stability, as will be discussed in Section 5.3. 



4.3. Spin problem for L > 0.3LEdd 
As mentioned in Section 1, some works on the black hole spin measure ment showed that 



the standard thin disk m odel is likely to be inaccurate for L > 0.3LEdd (e.g.. lMcClintock et al. 



20061 ; IStraub et al.l 1201 ll ) . One explanation is that the inner disk edge is still located at the 
innermost stable circular orbit (ISCO), but its emission is shaded by the outer disk. Thus, 
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the in ner disk radius obtained from the spectral fitting is not true. However, IWeng &: Zhang 



( 120 111 ) showed that the disks in black hole and neutron star X-ray binaries trace the same 
evolutionary pattern for L > 0.3L Edd . In addition, for the neutron star system XTE J 1701- 
462, the bo undary emissio n area maintains nearly constant despite the varying luminosity 



of the di sk (ILin et al.ll2009l Fig ure 17), which indicates that the neutron star's surface is not 
shaded. Weng fc ZhangJ ( 2011 ) therefore argued that the inner disk of the black hole system 



should not be shaded either due to the similar phenomenon. They suggested that the inner 
disk radius moves outward because of the increasing radiation pressure. 

In our opinion, from the energy advection, it is easy to understand that the standard 
disk model seems to be inaccurate above 0.3LEdd- As revealed by the lines in Figure 4, / a( j v 
is likely to be non-negligible (probably ~ 0.1) for m ~ 0.3 at r = 10r g . We would point out 
that, compared with the Paczyhski-Wiita potential, the Newtonian potential in the present 
work may magnify the viscous heating rate at small radii such as 10r g , thus the real / a d v at 
10r g may be smaller than the values showed in Figure 4. On the other hand, for the same 
m, since the viscous heating rate at a smaller radius such as r = 5r g will probably be larger 
than that at r = 10r g , so does the advection factor. We can therefore expect that, even 
for the Paczyhski-Wiita potential, the advection at the position close to the ISCO should 
be non-negligible for m ~ 0.3. Consequently, the standard thin disk model, based on the 
energy balance between the viscous heating and the radiative cooling with the advective 
cooling being ignored, may be inaccurate. 



4.4. Vertical height 

Figure 2 shows that p decreases significantly with decreasing 9, and Figure 4 implies 
that most of the accreted matter exists in the range n/A < 9 < 3n/A. In order to have a 
more clear view, we define an averaged dimensionless height as A9 = S/2rpo ; where the 
surface density £ takes the form: 

7T-0Q 

psin9 d9 . (20) 

Figure 5 shows the variation of / a d v with A9 (solid line). Even though the photosphere 
is close to the polar axis, the averaged height A9 is geometrically slim with 0.3 < A9 < 
0.6. Furthermore, the figure shows that / adv increases with increasing A9 or m, which 
agrees with the classic picture. For quantitative comparison, we plot the function / a d v = 
1.5tan 2 (Afl) (dashed line ) in Figure 5 due to the relationship / adv > (H/R) 2 introduced by 



Abramowicz et al.l (119951 ). which is equivalent to / adv ^ tan 2 (A#) here. It is seen that / a d v is 
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not well proportional to tan 2 (A#). In the range 0.6 < m < 1.1 or 0.3 < / a d v < 0.8, however, 
we may regard the formula / a d v = 1.5 tan 2 (A^) as a rough approximation. 

5. Two-dimensional solutions and convective stability 

5.1. Two-dimensional solutions 

In Section 4, we focus on the solutions at a typical radius r = 10r g . In this section 
we will study the disk solutions for various radii. Since the vertical solutions are based on 
the assumptions of partial derivatives in the radial direction (presented in Section 2.1), it is 
necessary to derive vertical solutions for various radii to check whether these assumptions are 
self-consistent. Following the example solution in Figure 2, we study the two-dimensional 
solutions for m = 0.6 in the range 6r g ^ r ^ 12r g and < 9 ^ it/2. Figure 6 shows the 
radial variations of c s and E (solid lines) for five polar angles, i.e., 9 = 90°, 75°, 60°, 45°, 
and 30°. For comparison, the radial profile of Vk, which is proportional to r -1 / 2 , is shown in 
Figure 6(a), and an example slope of oc r~ 5 / 2 is shown in Figure 6(b). The figure shows that 
c s and E behave roughly as oc r -1 / 2 and oc r~ 5 / 2 , respectively, which agrees with the original 
assumptions of radial derivatives. As mentioned in Section 2.1, once the radial derivatives 
of c s and E are given, the other ones can be inferred from Equations (l)-(7). Thus, we 
can expect that the radial derivatives of v r , v^, p, and F r in the two-dimensional solutions 
should also be in agreement with the assumptions. Our solutions are therefore likely to be 
self-consistent. 

In our calculation, the location of the photosphere does not vary much with the radius, 
i- e -> #o i$ 5° for various radii. As discussed in Section 4.1, the real position of the photosphere 
is likely to be located in the range 5° < 9 < 45°. We will make some comparison with 
simulations for the photosphere in Section 6. As shown in Figure 6, the derivatives of c s and 
E deviate a little for r — > 12r g and 9 = 90°. Such a divergence may be well understood from 
Figure 1, which shows that the solution for m = 0.6 and r — > 12r g is quite close to the lower 
solid line, which indicates that the gas pressure may not be negligible. In particular for the 
equatorial plane, the mass density has the maximal value there, thus the gas pressure may 
be most significant at this position. 

5.2. Solberg-H0iland conditions 

In this section, we will study the convective stability of the radiation pressure-supported 
disks. The well-known Solberg-H0iland conditions in cylindrical coordinates (R, <fi, z) take 
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the forms (e.g.. iTassoull 120001 ): 



i_az 2 

R 3 dR 



-T-r- - — — VP • VS > 



C P p 



dP 


(dl 2 dS 


dl 2 dS 


dz 


[dR dz 


dz dR 



> 



(21) 
(22) 



where I is the specific angular momentum per unit mass, P is the total pressure, Cp is the 
specific heat at constant pressure, and S is the entropy expressed as 



dS oc d In 



P 

7' 



(23) 



where 7 is the adiabatic index. 

The R and z components of the well-known Brunt- Vaisala frequency are written as 

P 



*r2 1 9P d , 

Nr = — —In 

R jp dR dR 

, T 2 1 9P d , 

^ = 7T7T ln 

7P a2; 



P 7 
P 7 



and the epicyclic frequency takes the form: 



K 



1_dP 

R 3 dR 



Thus, the first Solberg-H0iland condition, Equation (21), can be simplified as 



N% = N 2 r + N 2 + k 2 >0 



(24) 



where N e g is defined as an effective frequency. For accretion disks, there usually exists 
dP/dz < (as shown in Figure 8), so the second Solberg-H0iland condition, Equation (22), 
reduces to 

P\_dl^_d_ 
pV ~dzdR 

In numerical calculation, we adopt P = E/3 and 7 = 4/3 according to the radiation pressure- 
supported assumption. 



dl 2 d , 
AlS =dRd-z ln 



In 



P 

P 7 " 



> 



(25) 



5.3. Convective stability 



Based on the two-dimensional solutions for m = 0.6 in Section 5.1, we can obtain 
the variations of physical quantities in cylindrical coordinates and therefore investigate the 
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convective stability by Equations (24)-(25). We take the cylindrical radius R 
typical position to study the convective stability. 



10r g as a 



Figure 7 shows the ^-direction variations of k 2 , N r , N 2 , N 2 s , and Ai$, where the former 
four quantities are normalized by Q^, and the last one is normalized by v\. The positive 
values for both N 2 S and Ais indicate that the disk should be convectively s table. For the 



equato rial plane, the result of N 2 S > can be inferred by Equation (15) of iNarayan fc Yi 
(11994J ). which revealed that the disk will always be convectively stable for 7 = 4/3 at z = 0. 



For further understanding the convectively stable results for z > 0, we plot Figure 8 to 
show the z-direction variations of p, E, I, E/p^ 3 , and the advection factor / a a v . It is seen 
that p drops faster than E with increasing z, so the quantity E/p^ 3 increases with z. From 
Equation (23) we immediately have 

E \ 

(26) 



dS d 
dz dz 



,4/3 



> 



which is known as the Schwarzschild criterion for a constant angular velocity at a cylindrical 
surface (dQ/dz = 0), corresponding to the so-called barytropic flows where the pressure 
depends only on the density. As the dashed line in Figure 8 shows, the angular momentum 
I (or equivalently the angular velocity O) does not vary significantly with z. Thus, the 
convectively stable results are easy to understand from dS/dz > 0. 

In a similar study (vertical structure based on the local energy balance) of the general 
model for optically thick di sks, however, the disk was found to be convectively unstable (e.g., 
Sadowski et al.ll2009l . 120111 ). As interpreted by the three possible reasons in Section 4.2, the 
significant difference in the results between their works and ours is probably related to the 
different approach to describing the vertical structure. In addition, here we would mention 
two more details which may help to understand the difference in convective stability. First, 
the profiles of p and E in Figure 8 reveal that the absolute value of radial velocity \v r \ 
increases with z (\v r \ oc E/p inferr ed from Equations (7) and (10)). Compared with the 
uniform vr in lSadowski et all (120111 ). our increasing \v T \ with z may result in faster drop of p 



in the z direction (steep er slope of p) if we sim ply assume the mass supply to be comparable. 
Second, compared with ISadowski et al.1 (120111 ). the advection in our solutions is significantly 
stronger, which means that for the same rh thus comparable viscous heating rate, the vertical 
radiation flux F z will be less in our results. Thus, E may decrease slower in the z direction 
(flatter slope of E) due to the less F z and lower p (inferred from Equation (6)). As indicated 
by Equation (26), the flatter slope of E and the steeper slope of p will both make contribution 
to dS/dz > 0, and the disk is therefore likely to be convectively stable. 



We would stress that, our solutions are limited by the radiation pressure-supported case. 
Thus, the present results cannot directly show the convective stability of disks for either the 
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gas pressure-supported case or the case of comparable gas and radiation pressure. Actually, 
we have made some additional calculation for t hin disks in cylindrical coordinates to check 
the convective stability, following the method of Sadowski et al. (j2011 ) but without consid- 
ering advection. We found that the disk is convectively stable for the gas pressure-supported 
case, whereas the disk is convectiv ely unstable for t he ca se of radiation pressure being signif- 
icant. Thus, we would agree with ISadowski et al.l ( 120111 ) on the convectively unstable disks 
for moderate accretion rates such as 0.01 ~ 0.1M Edd , corresponding to significant radiation 
pressure and non-negligible gas pressure. Moreover, for m < 0.1, the disk will be geometri- 
cally thin, thus there is no difference between the assum ptions vp = and v z = 0, and / a a v 
is probably negligible. As a consequence, the solutions of 
accurate. 



Sadowski et al.l (120111 ) ought to be 



Furthermore, as revealed by the profiles of E and p in Figure 8, d In E/d In p is less 
than unity in the z direction. Following the argument in Section 4.1, the polytropic relation 
seems unsuitable either in the z direction. In addition, as shown in Figures 7 and 8, our 
example solution for rh = 0.6 at R = 10r g is terminated at z/R = 0.66. The reason is that 
the solutions in spherical coordinates is limited by r = 12r g (as shown by the horizontal 
dot-dashed line in Figure 1), which corresponds to z/R = 0.66 at R = 10r g in cylindrical 
coordinates. 



5.4. Analysis of convective stability for z <C R 

For the region close to the equatorial plane, we can make some analysis of the convective 
stability by the Taylor expansion method. Obviously, we have dS/dz = at z = from 
symmetric conditions. Thus, for z <C R, the value of dS/ dz can be estimated by the second- 
order derivative at z = 0: 

dS d 2 S , 

^«z w (,«R). (27) 

Based on Equations (7)-(12), we can eliminate v r , E, and Fq and therefore obtain a 
set of three equations for the three quantities t>^, c s , and p. By using the Taylor expansion 
method, together with the boundary condition of Equation (14), we derive the following 
three relationships for the second-order derivatives of v^, c s , and p (the ram-pressure term 
in Equation (8) is ignored): 

5c s d 2 c s d 2 v <t> 

— — + Va, — = , 

1 d 2 p | 2 d 2 c s _ v% 
P d~Q 2 c s d~Q 2 c 2 ' 
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3_d\ t _l_d^ji ^^Cs _J ( d 2 c s \ 

where 9 is defined as 9 = n/2 — 9, which is a small value in the analysis. 

In our solutions we have C ^ at z = 0, e.g., c^/v^ = 0.16 for m = 0.6 at R — 10r g . 
For the simple case with c 2 /i>| <C 1, the above three relationships provide 

d^nvtp _ h_v\ _ d 2 \nc s _ \ v\ _ <9 2 lnp ^ 

The second-order derivative of entropy at z = can therefore be derived by the following 
coordinate transformation: 

a's a 2 , / c, 2 \ i a 2 , / c j \ u /<j\ i /n>j a 

a? * a? ln (j^J = S 5 ^" 1 (y*J + 5* (y*J " F (4 - 2 J > ' (28) 

Thus, Equations (27)-(28) indicate dS/dz > for the region close to the equatorial plane. 
The disk in this region is therefore likely to be convectively stable. 



6. Summary and discussion 

In this paper, we have studied the vertical structure, energy advection, and convective 
stability of radiation pressure-supported disks in spherical coordinates. In the 9 direction, 
we replaced the pressure-density polytropic relation by the local energy balance per unit 
volume between the viscous heating and the advective cooling plus the radiative cooling, 
and obtained the distribution of physical quantities such as p, v r , v^, c s , E, and F e . The 
photosphere was found close to the polar axis and therefore the disk seems to be extremely 
thick. However, most of the accreted matter exists in a moderate range around the equato- 
rial plane such as n/A < 9 < 3n/A. We showed that the polytropic relation is unsuitable for 
describing the vertical structure of radiation pressure-supported disks. More importantly, we 
found that the energy advection is significant even for slightly sub-Eddington accretion disks, 
which is quite different from the previous result that the advection is of importance only for 
super-Eddington accretion disks. We argued that, the non- negligible advection may help to 
understand why the standard thin disk model is likely to be inaccurate for L > 0.3LEdd- In 
addition, we studied the two-dimensional solutions to check our basic assumptions of radial 
derivatives, which indicates our solutions to be self-consistent. Furthermore, we investi- 
gated the convective stability of the disks in cylindrical coordinates by the two-dimensional 
solutions derived in spherical coordinates. The disk solutions satisfy the Solberg-H0iland 
conditions, which reveals that the disk ought to be convectively stable. 
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To our knowledge, there are mainly two series of simulations on optically thick accre- 
tion flows. One is three-dimensional, radiation magne tohydrodynamic (RMHD) simulations 



by a shearing box, which focused on thin disks (e.g., iHirose et al.ll2006t iKrolik et all 120071 ; 
Hirose et al.ll2009al Jbl). The other is two-dimensional, either radiation hydrodynamic (RHD) 
or RMHD simu lations for global flows, which included thin disks and super-Eddington accre- 
tion flows (e.g.. lOhsuga et al.ll2005l l2Q09t lOhsuga fc Mineshigell201ll ). We would like to com- 
pare our numerical results with t he latter owing to th e global, radiation pressure-supported, 
and geometrically not thin case. lOhsuga et all ( 120051 ) studied RHD simulations in spherical 
coordinates with the a stress assumption for only the r<p com ponent, for extreme l y high 
accretion rates M = 300, 1000, and 3000L Edd /c 2 . Moreover, lOhsuga fc Mineshigel (boilh 
studied RMHD simulations in cylindrical coordinates for M ~ 100LEdd/c 2 . In our results, 
the position of the photosphere is quite close to the polar axis with < 5°. As discussed 
in Section 4.1, the real pos ition may exist in the range 5° < 0q < 45°. In the simulations, 
Ohsuga fc Mineshigel (120111) showed that the photosphere is around z/R = 2.4 correspond- 



ing to 6*o ~ 23°. lOhsuga et al.l (120051 ) did not mention the position of photosphere, but their 



Figure 4 for the density distribution reveals that the photosphere should exist at a certain 
#o significantly less than 45°. Thus, we would express that, for accretion rates around and 
beyond M^dd, the photosphere may exist far from the equatorial plane with 6$ < 45°. 

In addition, we discuss the possible link between our disk model and ultraluminous X-ray 
sources. As shown by the upper solid line in Figure 1, there exists a maximal accretion rate 
m miB varying with the radius. We argue that the possible upper limit of the accretion rate 
may help to understand why most ultraluminous X-r ay sources (ULXs) a re not in thermal 
dominant state. As mentioned in the review paper of iFeng fc Sorial (120111 ). there may exist 
three classes of black holes in ULXs: normal stellar mass black holes (~ 10M Q ), massive 
stellar black holes (< 100M o ), and intermediate mass black holes (10 2 — 10 4 M Q ). A massive 
stellar black hole with a moderate super-Eddington accretion rate seems to account for most 
sources up to luminosities ~ a few 10 40 erg s -1 . The slim disk model, which is the classic 
model for super-Eddington accretion disks, predicts a dominant thermal radiation from the 
disk. However, observations have shown that ULXs are not in the thermal dominant state 



excep t only a few sources such as M82 X-l (IFeng fc Kaaretll2010l ) and HLX-1 (IDavis et al. 
201 lh . In our understanding, the radiation of ULXs may be interpreted by an optically 
thick disk with fn < m max plus strong outflows. The disk will provide a thermal radiation, 
which is normally not dominant because of the moderate m max . On the other hand, the 
outflows m ay make contribution to th e non-thermal radiation by the bulk motion Comp- 
tonization ( jTitarchuk fc Zanniaslll998|) or through the jet of th e radiation-pressure driven 
and magnetically collimated outflow (lOhsuga h Mineshigel 1201 if ). 
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r/r g 

Fig. 1. — Solutions in the fn-r diagram. The parameter space is divided into three regions by 
two parallel solid lines. The middle region, denoted by "radiation pressure" , corresponds to 
the radiation pressure-supported disk, which is our main interest in this work. An example 
solution for rh = 0.6 at r = 10r g (filled circle) is shown in Figure 2. The solutions for 
various m at a typical radius r = 10r g (vertical dashed line) are focused on in Section 4. The 
two-dimensional solutions for m = 0.6 (horizontal dot-dashed line) are studied in Section 5. 
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Fig. 2. — Vertical structure of the disk for m = 0.6 at r = 10r g : (a) variations of p, \v r \, v^, 
and c s ; (b) variations of | i 7 ^ | / c£7 and t cs . 
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Fig. 3. — Vertical distribution of din E/ din p for m = 0.5 (dashed line), m = 0.6 (solid 
line), and m — 1 (dotted line). 




Fig. 4. — Variation of / adv with m (solid line) and 7^/4 (dashed line). 




Fig. 5. — Variation of / a d v with the dimensionless height A6 (solid line). For comparison, 
the function / a d v = 1.5 tan 2 (A^) is plotted (dashed line). The four typical accretion rates, 
m = 0.5, 0.6, 0.8, and 1, are denoted by filled circles. 
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Fig. 6. — Radial variations of c s and E (solid lines) for the polar angle 9 = 90°, 75°, 60°, 45°, 
and 30° for rh = 0.6. For comparison, the Keplerian velocity vk (oc r -1 / 2 ) and an example 
slope of oc r~ 5 / 2 are plotted by the dashed lines in (a) and (b), respectively. 
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z/R 

Fig. 7. — ^-direction variations of k 2 (dot-dashed line), (dashed line), N 2 (dotted line), 
N 2 G (thick solid line), and A/s (thin solid line) for m = 0.6 at a cylindrical radius R = 10r g . 
The quantities k 2 , 7V|, N 2 , and N 2 S are normalized by f^, and A/5 is normalized by v\. 
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z/R 

Fig. 8. — z-direction variations of (E / p 4 ^ 3 ) / (E / p^ 3 ) (thick solid line), l/l (dashed line), 
E/E (dotted line), p/po (dot-dashed line), and / a d v (thin solid line) for rh = 0.6 at R = 10r g , 
where the subscript "0" represents the quantities on the equatorial plane. 



